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ABSTRACT 

We investigate the stability of prograde versus retrograde planets in circular binary 
systems using numerical simulations. We show that retrograde planets are stable up 
to distances closer to the perturber than prograde planets. We develop an analytical 
model to compute the prograde and retrograde mean motion resonances' locations and 
separatrices. We show that instability is due to single resonance forcing, or caused 
by nearby resonances' overlap. We validate our results regarding the role of single 
resonances and resonances' overlap on orbit stability, by computing surfaces of section 
of the CR3BP. We conclude that the observed enhanced stability of retrograde planets 
with respect to prograde planets is due to essential differences between the phase-space 
topology of retrograde versus prograde resonances (at p/q mean motion ratio, prograde 
resonance is of order p — q while retrograde resonance is of order p + q). 



1 INTRODUCTION 



1 The stability of coplanar prograde planet orbits in 

■ binary systems h a s bee n investisated n umeri cally by 
; iHolman fc WiegertI (|l999l '). iMudrvk fc Wul (|2006l ) showed 

that instability in eccentric binaries is due to overlap of sub- 
resonances associated with certain mean motion ratios p/q. 
These sub-resonances are split due to the precession rate in- 

■ duced by the secondary star, hence overlap of sub-resonances 
for a given mean motion ratio p/g extends over a wide region 

, and expla ins the instability regio ns in eccentric binaries. Ad- 
' ditionaly, iMudrvk fc Wul (|2006l ) suggest that the cause for 
] instability in circular binaries is overlap of sub-resonances 
associated with the 3/1 mean motion ratio. However, this 
cannot work for circular binaries since in this case there is 
only one resonant angle. 

The circular restricted 3-body problem (CR3BP) is 
the simplest theoretical tool to understand planet stability 
within a binary system. The existence of an integral of the 
motion (Jacobi constant) reduces the number of variables of 
the coplanar problem from 4 to 3. Therefore, the phase space 
topology can be investigated by using surfaces of section for 
any given mass ratio fi. The Jacobi constant and associated 
zero velocity curves (ZVC) can impose bounds on the test 
particle's motion. In particular, it is useful to compare the 
Jacobi constant with the values at the coUinear Lagrange 
points (Z/i, L2 and L3). When the Jacobi constant exceeds 
the value at Li , the test particle must remain in orbit around 
eit her primary or s econdary stars (concept of Hill stabil- 
ity, [si^eheii iHiO) ) . When the Jacobi constant is smaller 
than the value at Li, the test particle can orbit both stars 
and will eventually collide with one o f them, although these 
capture episodes can be lo ng-lived l|Winter fc Vieira Netd 
I2OOII : lAstakhov et all 120031 '). When the Jacobi constant is 
smaller than the values at L2 or L3, the test particle can es- 



cape through these points. However, its is well known that 
the se are necessary but not sufficient conditions for instabil- 
ity (|Szebehelvlll980l). 
lEberle et al. 

I diooi) investigated stability of prograde 
planet orbits within circul ar binary systems, ba sed on the 
Jacobi constant criterion. iQuarles et all l|201lf ) validated 
these results by computing the maximum Lyapunov expo- 
nent which is a mea sure of chaos and associated instability. 
iQuarles et"ai](|201ll ) show that if the ZVC opens at L3 then 
the orbit is unstable but when the ZVC opens at LI or L2, 
the orbit is not necessarily unstable. The interpretation of 
these results is not obvious and it may depend on the par- 
ticular choice of initial conditions. 

IChirikwl (|l960l ) established the resonance overlap crite- 
rion to explain the on set of chaotic m otion in Hamiltonian 
deterministic systems. IWisdomI (|l980l ) obtained a resonance 
overlap criter ion for th e onset of chaos in the CR3BP valid 
when < 1. IWisdomI (jl980,) showed that first order mean 
motion resonances with the secondary overlap in a region of 
width ~ fi^^'^ . Orbits in this region exhibit chaotic diffusion 
of eccentricity and semi-major axis until escape or collision 
occurs. When fi <^ 1 individual resonances cannot increase 
the eccentricity or semi-major axis up to escape values. How- 
ever, in binary star systems ~ 1 thus the effect of single 
resonances is not necessarily negligible. 

Recently, it was possible to detect the Rossiter- 
MacLaughin effec t on transiting extra-solar planets 
(|Triaud et al.ll2010D . This effect allows to measure the ori- 
entation of the planet' orbit with respect to the parent 
star's equator. Contrary to what happens in the Solar 
System, extra-solar planets' orbits can be misaligned with 
the host star's equator with angles that range from 0° to 
180°. Several mechanism have been proposed to explain 
these misaligned planets. These include classic Kozai os- 
cillations due to a nearby star with subsequent tidal drift 
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l|Correia et alj I2OIII I. secular interaction with a compan- where 



ion brown dwarf or giant planet followed by tidal drift 
(jNaoz et al."201l''). secular chaos and tides in multiple planet 
systems (|Wu fc Lithwick 2011,), orbit inte raction in multipl e 
planet systems with a companion sta r (iKaib et al. 201l[). 
or pl anet-planet scattering and tides l|Beauge fc Nesvornvl 

|201]J). ^ ^ 

'^Ga von fc Boi^ (2008) used the MEGNO chaos indica- 
tor (,Cincotta fc Simo 2000 ) to show that retrograde reso- 
nance in 2 planet systems is more stable than the equiva- 
lent prograde resonance. Therefore, they suggest that ret- 
rograde resonance could explain the radial velocity data of 
extra-solar systems where prog rade resonance is unstable 
l|Gavon-Markt fc Bold [20091 '). In lGavon et~aLl (|2009l '). an ex- 
pansion of the Hamiltonian for retrograde resonance in 2 
planet syste ms is deve l oped. However, the numerical explo- 
ration in G avon et al.l (|2009| ) is limited to a small set of ini- 
tial conditions which could explain why they do not conclude 
on the essential differences between prograde and retrograde 
resonance. 

Planets in retrograde orbits within a binary system are a 
theoretical possibility although none was confirmed to date. 
A planet on a retrograde orbit has been suggested as an ex- 
planation for a per iodic signal of 416 d in i^-Octantis A ra- 
dial velocity curve l|Ramm et al.|[2009h . The star I'-Octantis 
A orbits its companion (z^-Octantis B) on a 2.9 yrs orbit. 
Such a tight binary orbit implies that a prograde planet at 
416 d is unstable but a retrograde planet c ould be stable 
at least up to lO** yrs (lEberle fc Cuntj|2010l ). Nevertheless, 
there are alternative hypothesis that claim that I'-Octantis 
A radial velocity could be explained, wit hout the need of a 
plane t, if !/-Octantis B was a double star (jMorais fc Correial 
I2OI2D . 

The purpose of this article is to investigate stability 
of coplanar prograde and retrograde planet orbits in circu- 
lar binary systems. Contrary to previous works, we will not 
only perform simulations but will also provide theoretical 
explanations for the onset of instability based on the effect 
of single resonances or due to resonance overlap. 



2 EXPANSION OF THE DISTURBING 
FUNCTION IN CR3BP 

We consider the planar CR3BP composed of a test parti- 
cle orbiting a primary mo, and perturbed by a secondary 
7712. The primary and secondary h ave a circular orbit with 
frequency 712 — \/ G{mo + m2)la\ and radius a2. Since we 
want to model the perturbation from the secondary we write 
the equation of motion in the frame with origin at the pri- 
mary 



ri = -V (C/o + U) 



(1) 



where l7o = Gmo/ri, G is the gravitational constant, and 
U is the disturbing potential due to the perturber m2. 

When U — Q (i.e. 7712 = ) the solution to Eq. (1) 
is a Keple rian elliptical orbit with mean motion 77i = 
^J G mo /a\, semi- major axis ai, eccentricity ei, longitude 
of pericenter vji, and true anomaly /i. 

The disturbing potential is 



U = Gm,2 



/I a r\ 
\A 02 ai 



cos S 



A = ||ri — r2|| = J + a| — 2 ri a2 cos S 



(3) 



a = ai/a2 < 1, and S is the angle between ri and r2. 

In the prograde case the primary-secondary and test 
particle orbit in the same direction. In the retrograde case 
the test particle orbits in the opposite direction of the 
primary-secondary. The primary-secondary relative posi- 
tion vector is r2 = 02 (cos(A2), sin(A2)).The test particle 
{nil ~ 0) position vector with respect to the primary is 
n = n (cos(/i-|-7^i),sin(/i+ii7i)) with n = ai (l-e?)/(l + 
ei cos /i). Hence, 



cos 5* — cos(/i + ZUi — A2) 



(4) 



where A2 ~ ±712 t, and the ± sign applies to the prograde 
or retrograde cases, respectively. 

The 1st and 2nd terms in Eq. (2) are known as di- 
rect and indirect parts, respectively. The disturbing po- 
tential (Eq. (2) can be expressed in the orbital elements 
(ai, ei, /i, cc7i) by using Laplace coefficients (literal expan- 
sion). The direct part is written as a Taylor series in e = 
(ri/m - 1) i.e. 



^= 1 + 



1 « « rfM 1 



(5) 



with 

1 

P 



0.2 



-{1 + a ~ 2a cosS) 



-1/2 



i-5]i^/,(a)cos(,S) 



a2 ^ 2 

i 



(6) 



where IP-^^^^a) is a Laplace coefficient. 
Since 

cos(j5') = cos(j(./'i +t^7i A2)) , (7) 

using elliptic expansions for ri/ai, cos/i and sin/i, we ob- 
tain for any give n ,7, the direct and indi rect parts of Eq. ((2]). 
This is done in lElfis fc Muirwl (|2000[ ) for prograde reso- 
nances. In the planar CR3BP, we only cons ider those terms 
in the expansion from lEllis fc Murravl (|2000l ') that depend on 
ei . These terms consist of cosines of angles which are combi- 
nations of the mean longitudes A2 = 712 t, Ai = rii (t— r)-f lui 
(where r is the time of passage at pericenter), and the lon- 
gitude of pericenter roi. From the discussion above we con- 
clude that for retrograde resonances the terms are exactly 
the same, although we must replace A2 = —712 t- 

B y inspecting the expa nsion of the disturbing poten- 
tial in lEllis fc Murravl (|2000l ) we see that at 1st order in ei, 
terms of the type (,?' — 1) Ai — j A2 -I- vJi appear (4D1.1 in 
lEUis fc Murravl (|2000l ')'). If j > 2, these terms correspond to 
^ ~ 1) prograde resonance since Ai = ni and A2 — 712 
thus the time variation of the angle is {j — 1) 7ii — j 712 ~ 0. 
In the retrograde case, A2 = —712 hence the previous terms 
are non-resonant. At 2nd order in ei ter ms of the type 
7 - 2 ) Al - j X2 + 2zui appear (4D2.1 in lEllis fc Murravl 



(2) 



200y)) which correspond to a — 2) prograde resonance 
^ 3) or the 1/1 retrograde resonance when j = 1. At 
3rd order in ei t erms of the typ e (j — 3) Ai — j A2 + 3 uji 
appear (4D3.1 in lEUis fc Murravl (|2000l )') which correspond 
to a j/{j — 3) prograde resonance {j > 4) or the 2/1 retro- 
grade resonance when j = 2. At 4th order in ei terms of the 
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type ( .7-4) A1-.7 A2+4roi appear (4D4.1 in lEUis fc Murray! 
(|200CI )) which correspond to a — 4) prograde resonance 
> 5) or the 3/1 retrograde resonance when j = 3. It 
can be shown that at 5th order in ei terms of the type 
{j — 5) Ai —j A2 +5 ci7i appear which correspond to a j/{j — 5) 
prograde resonance {j > 6) or to the 3/2 retrograde reso- 
nance when j — 3 and the 4/1 retrograde resonance when 
j = 4. Therefore, we see that p/q prograde resonances are 
of order p — q while p/q retrograde resonances are of order 
p + q. The only resonant terms in the indirect part of the dis- 
turbing function (C R3BP) corres pond to the 1/1 prograde 
resonance (4E0.1 in lEllis fc Murr ay (200( |)) and to t he 1/1 
retrograde resonance (4E2.2 in lElhs fc Murray! (HooO)). The 
secular term is obtained by averaging over the mean longi- 
tudes Ai and A2, hence i t is the same in the prograde or 
retrograde case (4D0.1 in !Enis fc Murravl ^^)). 



3 ANALYTIC MODEL FOR MEAN MOTION 
RESONANCE IN CR3BP 

Here, we briefly review the analytic model for 1st, 2nd 
and 3 rd order prograde resonance from !Murray fc Dermott] 
1! 19991 ). In Appendix A we derive in detail this Hamilto- 
nian model in the framework of the CR3BP. We extend the 
model to retrograde resonance of lowest (3rd) order. We ex- 
plain how we can use the model to obtain resonance widths 
for initially circular orbits. 



3.1 Hamiltonian model for prograde/retrograde 
resonance 

In the CR3BP, j/ {j — k) prograde resonance is of order k 
while j/{k — j) retrograde resonanc^ is of order k. The 
resonant angle is 



( j - fc) Ai - j A2 + fc zui 



(8) 



Here, we will summarize the results for prograde ~ 
k) resonances of 1st, 2nd and 3rd order (fc = 1,2,3) and 
we will extend these results to the retrograde 2/1 resonance 
which is of 3rd (lowest) order (j = 2, fc = 3). The reso- 
nant Hamiltonian (Eq. IA18|) depends on a single parameter 
(F,a. [Kw\ 



S — A[(j — k)nl T j"-2 + fc'd7i]/fc , 
where n1 — ni + Xl, 

/24-fc j*=-s/3(j _ fc)'=— */^fc'*~2'=\ ^ 

/i2/d(Q)2 

and from Eqs. (IA7IA8|l with ei < 1 

uji « 2 afs(Q-)ni 

mo 

i. ^ ■ . 



(9) 

(10) 

(11) 
(12) 



with values of fs{ce) at resonant a — [j/\j — k\] shown 
in Table 1. 



^ We will use the notation j/ (k—j) retrograde resonance or j/ (j — 
k) resonance: e.g. 2/1 retrograde resonance or 2/-1 resonance 
(j = 2, fc = 3). 



resonance 


a 


a fs (a) 


a/d(") 


4/1 


0.39685 


0.032355 


-0.09698 


3/1 


0.48075 


0.068381 


-1-0.28785 


5/2 


0.54288 


0.11600 


-0.61503 


2/1 


0.62996 


0.24419 


-0.74996 


2/-1 


0.62996 


0.24419 


-0.25304 


5/3 


0.71138 


0.51566 


-f2. 32892 


3/2 


0.76314 


0.87975 


-1.54553 



Table 1. Values of secular an resonant functions at resonant a ■ 

ij/\j-k\)-y^,k = 1,2,3. 



The Hamiltonian (Eq. IA18|I is expressed in cartesian 
canonical variables 



X = R cos{6/k) 
y = R sm{e/k) 

where the scaling factor is (Eq. IA20p 



R 



3(-l)' (3 



fc)4/3j2/3 



.fd{a)fi 



fc2 



er 



(13) 
(14) 



(15) 

fcir^'''* shown 



with values of fd{ct) at resonant a — 
in Table 1. 

Obviously, this analytic model is valid only for small 
perturbation i.e. (Eqs. IA1IA2IA3) 



Ures 1 "12 r , ^ 

ho 2 mo 



< 1 



Usee 1 "12 r r \ 2 . 

IT = -— Q.fs(a)ei < 1 
Uo I mo 



(16) 
(17) 



In Table 1 we show the values of a fd{a) and a fs{a) at 
resonant value a — — k\]~^^^. These provide a measure 
of the analytic model validity. In particular, the resonance 
location {S = 0, Eq. [9]) depends on the secular term {vji, 
Eg. Ilip . hence it is only accurate when Eq. (|17p is verified. 
From Table. 1 we see that the secular term ( Eq. (I17p ) at 
the 4/1 resonance is about 0.47, 0.28 and 0.13 times that 
at the 3/1, 5/2 and 2/1 (or 2/-1) resonances, respectively. 
Therefore, the secular term (Eq. (I17|l ) is approximately the 
same when fi = 0.2, /i = 0.09, /x = 0.07 and pi, = 0.03 at the 
4/1, 3/1, 5/2 and 2/1 (or 2/-1) resonances, respectively. 



3.2 First, second and third order resonances 

The Hamiltonian (Eq. lAlSP when fc = 1 is 



(18) 



When (5 < —3 there is a single stable equilibrium point 
(Fig. (iK). At (5 = —3 an unstable equilibrium point appears 
which bifurcates into a stable/unstable pair visible when 
(5 < -3 (Fig. [lb). 

The Hamiltonian (Eq. IA18|) when fc = 2 is 

ii2 = \ {x' + y') + l {x' + y^f + 2 (x^ - y^) . (19) 

At (5 = 4 the origin becomes an unstable point and 2 stable 
points appear at = ±7r/2 which move away from the origin 
as 5 increases (see 5 = (Fig. (2^) and 5 = — 4 (Fig. HJj)). 
At 5 = —4 the origin becomes again a stable point (Fig. [2]d) 
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Figure 1. Curves of constant Hi: 6 = (a) and S = —3.78 (b). 
The separatrix intersects the origin at (5 = —3.78. 



and 2 unstable points appear at = 0, vr that are visible 
when S < -4 (Fig.[2t). 

The Hamiltonian (Eq. IXTS)) when fc = 3 ifl 

H,= ^-{x'+y') + l{x' + yY -2xix^~3y^) . (20) 

At 5 = 9 the origin is a stable point and 3 pairs of sta- 
ble/unstable points appear at = 0,±2 7r/3. The 3 stable 
points move away from the origin while the 3 unstable points 
move towards the origin (Fig. [3^), until they coincide with 
it at J = (Fig. [3}3). At 5 = (exact resonance) the ori- 
gin bifurcates into a stable point and 3 unstable points at 
(j) = ±7r/3,7r. These unstable points move away from the 
origin as 5 decreases (Fig. [3];). 



6 n 




^ The diagrams with curves of constant H3 in lMurrav fc Dermotti 
lll999l l should be rotated by tt/3. 



Figure 2. Curves of constant H2'- 5 = (a), 5 = —4 (b) and 
5 = —6 (c). The separatrix intersects the origin between 5 = 
and 5 = —4. 
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(a) 



(b) 



(c) 

Figure 3. Curves of constant H3: <5 = 8 (a), 5 = (b) and 
5 = —8 (c). The separatrix intersects tlie origin only at (5 = 
(exact resonance). 



3.3 Computing the resonances' separatrices 

In order to comput e the resonance s' widths we follow the 
method described in IWisdomI l|l980l) for 1st order mean mo- 
tion resonances. This method was developed specifically for 
initially circular orbits and does not rely on the pendulum 
approximation which is more appropriate near the resonance 
center at ei 7^ 0. Wisdom's method consists in measuring 
the variation in the parameter, i.e. AS, between exact res- 
onance {S = 0) and the last value at which the separatrix 
intersects the origin (i.e. when the orbit with ei = is at 
the separatrix). From Eq. ([9}, we have 



\AS\ =A- 



■Am 



(21) 



and using Kepler's 3rd law we obtain the resonance half- 
width 



(22) 



Aai _2k \ A5\ 
ai 3j A 

For 1st order resonances the separatix intersects the 
origin at 5 ~ —4 (Fig. [T][b)) while for 2nd order resonances 
the separatrix intersects the origin from 5 = (Fig. [2]^ a)) to 
5 = —4 (Fig. [2l^b)). Hence, we take a range \ AS\ « 4 and ap- 
ply Eq.[22]to obtain the resonances' half- widths. We checked 
that the 1st order resonances' width s are in agreem ent with 
Wisdom's approximate expressions (|Wisdomlll980l ). 

For 3rd order resonances the separatrix only intersects 
the origin at 5 = (Fig. EJb)). Hence, AS = i.e. 3rd or- 
der resonances have zero width. This predicted zero width 
at ei = is a known feature of the analytic models, e.g. 
it also occurs when using the pendulum appr oximation 
iMurrav fc DermottI (ligogP : iMudrvk fc Wul (|2006D . 

The exact resonance location is obtained by solving S = 
(Eq.|9} for a = 01/02. In Sect. 5 we will see that, at small 
to moderate fi values, the resonance's widths / locations are 
in reasonably good agreement with the numerical results 
obtained by the method of surfaces of section. 



4 INITIAL CONDITIONS AND ZERO 
VELOCITY CURVES 

We chose a binary system with masses mo = Mq (primary) 
and m2 < Mq (secondary), inter-binary distance 02 = 1 AU, 
and mass ratio fi = m2/(mo+m2). We chose units such that 
G {mo + m2) = 1 which implies a binary period T2 = 2tt yrs. 
The initial orbital elements with respect to the primary were 
semi- major axis ai, eccentricity ei = 0, mean longitude 
Ai = 0, inclination / = (prograde orbits) or I = tt (ret- 
rograde orbits). Hence, the test particle was always started 
between the primary and secondary, orbiting in the same 
direction (prograde orbit) or in the opposite direction (ret- 
rograde orbit). 

The planet's orbit can remain bounded to the primary 
(stable orbit), or it can become unstable. Unstable orbits col- 
lide with either primary or secondary, or escape from the sys- 
tem. We assume collision with the primary if ro < 0.005 AU 
(i.e. the test particle gets within one solar radius of mo), 
collision with the secondary if ro < 0.005 m2/mo AU, and 
escape from the system if r > 3 AU. In practice, temporary 
capture in chaotic orbits around the secondary is possible 
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but all these capture episodes end either by collision with 
the stars or escape from the system. 

The CR3BP describes the motion of a test parti- 
cle in the frame c o-rota tinp; with the binary (see e.g. 
iMurrav fc DermottI (|l999l )l. In our problem the test par- 
ticle moves in the same plane as the binary and the position 
vector with respect to the binary's center of mass has coor- 
dinates {x,y). The CR3BP has a n integral of motion known 
as Jacobi constant l|Murrav fc Dermott 1999') 



r2 



■ 2 .2 

X -y 



where 



2 
2 

r2 



+ y 



(23) 

(24) 
(25) 



Due to our choice of initial conditions we have y{0) = 
and i(0) = 0, hence we can visualize the orbits using 
surfaces of section i.e. plotting {x, x) when y — and y x 
y{0) > 0. Since the test particle has ei = at t = then 



C = x{0f + 2 



1-M 



+ 



x(0) + M l-Ai-x(O) 



with 

x{0) 



m = ± 



1-M 

ai 



ai 



^)(0)^ (26) 

(27) 
(28) 



where the ± sign in y{0) applies to prograde and retrograde 
orbits, respectively. 

The zero velocity curves (ZVC) are obtained by solv- 
ing Eq. ^ with «2 = + = 0. These ZVC provide 
boundaries on the test particle's motion since it can only 
occur in the region with > 0. In particular, the ZVC at 
the coUinear Lagrange point Li is the limit curve for mo- 
tion solely around primary or secondary. The ZVC at the 
coUinear Lagrange points L2 and L3 are the limit curves 
that prevent escape from L2 or L3, respectively. Therefore, 
comparing the test particle's Jacobi constant with the values 
at Li, L2 and L3 provide us important information regard- 
ing stability. If C > Ci the test particle must remain in orbit 
around the primary. If C < Ci collision with secondary or 
primary stars is possible. If C < C2 or C < C3 escapes are 
possible through L2 or L3, respectively. 

Lagrange points have x = y = 0. The coUinear Lagrange 
points have y = while th e coordinate x can be ob tained 
as series expansions m fj. (|Murrav fc DermottI [l999l ) . The 
Jacobi constant at Li, L2, -L3, including terms up to 2nd 
order in fj. is 

^ „ , „4/3 2/3 10 

Oi 



3 + 3*/V'^' 



■ /i + - 3 ' 11 



2/3 4/3 



52 1/3 5/3 , 62 2 



C2 



C3 



14 1 
— P. + - 



32/3 ^4/3 



56 1/3 5/3 , 98 2 
"81^^ +81^^ 
01 1 2 



(29) 

(30) 
(31) 



In the next section we will plot the initial conditions that 
have C — C'l, C = C2 and C — C3 (where C is given by 



Eq. p6|) '). We will see that, as expected, these curves sepa- 
rate the regions of different end states for the test particle. 



5 NUMERICAL STABILITY STUDY 

We constructed grids of initial conditions in the plane {a, jj.) 
with an step A/i = 0.002 and Aa = 0.005 AU. Each point 
in the grid was then numerically integrated over 50 ~ kyr 
(around 12000 binary periods, depending on fi) using a 
Burlisch-Stoer based N-body code (precision better than 
10~^^) using astrocentric osculating variables. During the 
integrations we computed the averaged MEGNO chaos in- 
dicator (Y) (MEGNO is the acronym of Mean Expo nen- 
tial Growth of Nearby Orbits) (|Cincotta fc Simdll2000l ). We 
show these MEGNO maps in Figgfa) and Fig. Oa). 

The MEGNO chaos maps use a threshold that is cho- 
sen in order to avoid excluding stable orbits that did not 
converge to their theoretical value or those orbits that are 
weakly chaotic. Thus the color scale shows "stable" orbits 
in blue up to (Y) « 2.0 (a particular choice based on inte- 
gration of individual orbits for very long times and due to 
the characteristics of this system). 

MEGNO is a fast chaos indicator that allows to distin- 
guish rapidly between regular and chaotic orbits. Within the 
integration time, all the orbits identified as unstable (red) 
either collide with primary or secondary, or escape from 
the system as we can see in Fig. gjbfcc) and Fig. [Sfbfcc) . 
The thin region colored with light-blue in the transition be- 
tween chaos/regularity is typically unstable within ~ 10000 
to 15000 binary periods. We integrated the grids with less 
resolution for 2.5 x 10"" binary periods and no additional 
signs of instability were observed. 



5.1 Prograde case 

In Fig. |4] we show, for prograde orbits, the maps with: (a) 
MEGNO chaos indicator; (b) times of disruption of 3-body 
system; (c) planet end states (stable, collision or escape). In 
Fig. [5] we show a zoom of Fig. HJa). 

The separatrices of 1st order mean motion reso- 
nances (2/1 and 3/2) are shown as white dashed lines in 
Figs. |4ja)fc(b) and Fig. [5] The separatrices of 2nd order 
mean motion resonances (3/1 and 5/3) are shown as white 
solid lines in Figs. |3Ia)fc(b) and Fig. [S] These separatrices 
are obtained with Eq. (|22|) . The locations of 3rd order mean 
motion resonances (4/1 and 5/2), obtained by solving 5 — 
(Eq.[5| for a, are shown as solid grey lines in Figs. |3fa)fc(b) 
and Fig. [5] 

The initial conditions that have C — Ci and C — C2 
are shown as black solid lines in Fig.UJc). The test particle's 
end states depend on the the Jacobi constant. Orbits with 
C < Ci can escape through Li while orbits with C < C2 
can escape through L2. Since C > C3 escape through L3 is 
not possible. Moreover, the region near a = 1 has C > Ci 
hence escape is not possible. However, initial conditions in 
this region correspond to orbits in a collision route with the 
secondary at t = 0, thus they are unstable. 

From the perturbation theory, we predict that oscilla- 
tions in ei at the 2/1 resonance are large enough for collision 
with the secondary when /x>0.03, since an initially circular 
orbit at exact resonance reaches ei such that a(l + ei) ~ 1 
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Figure 4. Dynamical analysis for prograde orbits, a) Stabil- 
ity map in the (a, fi) space. Blue colors indicate stable orbits 
(< Y >< 2.0) while red colors indicate highly unstable orbits, b) 
Disruption times, c) End state of the planet. The separatrices are 
shown in (a)&(b) as dashed white lines (1st order resonances), 
solid white lines (2nd order resonances) and solid grey lines (3rd 
order resonances). The black solid curves in (c) indicate the initial 
conditions with C = Ci and C = C2- The white squares at sta- 
bility/instability transition zone indicate regions where we used 
the method of surfaces of section. 



Figure 5. Zoom of Fig. Hla). The separatrices are shown as 
dashed white lines (1st order resonances), solid white lines (2nd 
order resonances) and solid grey lines (3rd order resonances). The 
black solid curve a = 1 - 1.33 u^/'' is the 1st order resonance s 
overlap criterion l lWisdomlllQSOl ). 



(where ei is obtained from Eq.[T5]with R = 2, cf. Fig.[TJa)). 
This is in agreement with Figs. |4]&[5l 

We can see (Figs. |4] & [5]) that, at small to moderate 
fi values, the border of the unstable region seems to coin- 
cide with resonances' locations, namely the 4/1 resonance 
at a ft; 0.4 and the 3/1 resonance at a ft; 0.5. We can also 
identify chaotic regions that seem to be associated with res- 
onances' overlap (namely, between resonances 5/2 and 2/1, 
2/1 and 5/3, and 5/3 and 3/2). From Fig. [5] we see that 
the border of the unstable region when a > 0.7 approxi- 
mately coincides with the 1st order mean motion resonances' 
overlap c r iterio n (black curve in Fig. [5]) in agreement with 
IWisdonj (|l980l ). However, the perturbation method from 
which we obtained the resonance's locations and widths is 
certainly not valid when jj, is large or when a ~ 1. Therefore, 
we will plot surfaces of section, as described in Sect. 4, for 
initial conditions near the stability border (white squares in 
Fig. |4]). The integration time for the surfaces of section is 
10^ binary periods. 

Comparing the initial conditions in Fig. [4] (white 
squares) with their follow up orbits in Fig. [S] we see that per- 
turbation theory seems to be valid up to /i ft; 0.2 at the 4/1 
resonance. Following the discussion at the end of Sect. 3.1, 
we can extrapolate validity limits of ft; 0.09, fi ft; 0.07 and 
H ft; 0.03 at the 3/1, 5/2 and 2/1 resonances, respectively. 

In Fig. [HI a) instability is associated with the 4/1 res- 
onance separatrix {5 < case; cf. Fig. [Sj^c)). In Fig. [6jb) 
instability is also associated with the 4/1 resonance separa- 
trix (9 > 5 > case; cf. Fig. El^a)). In Fig. ^c) we identify 
a high order (k=22) resonance. In Figs. [6^d)&(e) instabil- 
ity is associated with the 3/1 resonance separatrix {S < —4 
case; cf. Fig. EI^c)). In Fig. [6l^f) we show stable orbits as- 
sociated with the 2/1 resonance (0.547 < a < 0.6). When 
a < 0.547 there is chaos due to overlap with 5/2 resonance. 
When Q > 0.6 eccentricity forcing at 2/1 resonance is large 
enough for collision with secondary, as seen above. 

We conclude that instability for prograde orbits is either 
due to single resonance forcing or resonance overlap. Regard- 
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(e) X (f) 

Figure 6. Surfaces of section near stability boundary, (a) /i = 0.165: Orbit witli a = 0.41 (red) shows inner circulation in 4/1 resonance. 
Orbits with a = 0.413 (magenta) and a = 0.414 (blue) are chaotic due to 4/1 resonance separatrix crossing. The blue orbit collides 
with the primary after 12740 binary periods, (b) fi = 0.25: Orbits with a = 0.45 (red), a = 0.46 (magenta) and a = 0.461 (blue) are 
in the vicinity of the 4/1 resonance separatrix. The blue orbit collides with the secondary after 8710 binary periods, (c) fi = 0.3: Orbit 
with a = 0.47 (red) is librating in high order (22) resonance and nearby orbit with a = 0.474 (blue) is chaotic. The blue orbit collides 
with the secondary after 1615 binary periods, (d) fi = 0.1: Orbits with a = 0.477 (red), a = 0.49 (magenta) and a = 0.497 (blue) are 
in vicinity of 3/1 resonance separatrix. The blue orbit collides with the secondary after 6800 binary periods, (e) /i = 0.051: Orbits with 
a = 0.51 (red) and a = 0.513 (blue) are in vicinity of 3/1 resonance separatrix. The blue orbit collides with the secondary after 25100 
binary periods, (f) fi = 0.051: Orbits with a = 0.547 (red), a = 0.59 (magenta) and a = 0.60 (blue). Stable orbits have 0.547 < a < 0.6. 
If Q < 0.547 there is overlap of 5/2 and 2/1 resonances. If a > 0.6 eccentricity forcing at 2/1 resonance is large enough for collision with 
secondary. 



ing the latter mechanism, we infer resonance overlap from 
the observation of a main resonance's chaotic separatrix. 
We know from Chrikov's criterion that widespread chaos in 
Hamiltonian systems is caused by resonances overlapping. 
However, we cannot always identify the resonance(s) that 



overlap with the main resonance. These are likely to be high 
order resonances which are difficult to identify in the sur- 
faces of section, in particular in the chaotic regions where 
the resonances' overlap. 

In Fig. [7] we present a case where we identify the over- 
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Figure 7. Proximity of 5/2 and 2/1 resonance overlap when /i = 
0.051 and a = 0.547. Evolution of 2/1 resonance angle (a) and 
5/2 resonance angle (b). 

lapping resonances. We show the evolution of the 2/1 and 
5/2 resonant angles when ^ — 0.051 and a = 0.547 (red 
orbit in the surface of section Fig. [6l^f)). The 5/2 resonant 
angle alternates between libration and circulation since it is 
at the separatrix. When a < 0.547 both separatrices over- 
lap and the orbits are unstable. These resonant angles are 
obtained from osculating elements with respect to the pri- 
mary. When the perturbation from the secondary is small 
the osculating elements are approximately Keplerian in the 
short term. Here, we can see that the resonant angles exhibit 
short term oscillations which indicates that the assumption 
of Keplerian osculating elements is not very good, despite 
the moderate mass ratio {fi — 0.051). Therefore, the oscu- 
lating elements cannot be used to identify the resonances at 
larger mass ratio jj, or when a ~ 1. The method of surfaces 
of section is always valid thus it is very useful to identify the 
main resonances and their chaotic separatrices. 



orbiting the secondary at t = 0, and is in agreement with 
the Jacobi constant criterion (C > Ci). 

In Fig. [S] we see that chaos and instability occurs at 
large values of the mass ratio fi or when a ~ 1. However, 
from the discussion at the end of Sect. 3.1 and the results 
in Sect. 5.1, we conclude that perturbation theory is valid 
only up to ~ 0.03 at the 2/-1 resonance. This threshold is 
well below the instability region which at the 2/-1 resonance 
occurs only when fi > 0.15 (Fig.fS]). Therefore, perturbation 
theory cannot be used and instead we will plot surfaces of 
section, as described in Sect. 4, for initial conditions near the 
stability border (white squares in Fig. [S)). The integration 
time for the surfaces of section is 10^ binary periods. 

In Figs. [9ja)&(b) instability is associated with the 2/- 
1 resonance separatrix (5 < case; cf. Fig. Efc)). As we 
noted above, in the retrograde case perturbation theory can- 
not be used to explain the instability border. In particular, 
the theoretical 2/-1 resonance location is over-displaced to 
the left as fi increases. This is mostly due to Eq. pip over- 
estimating the precession rate. Taking this into account we 
obtained a "corrected' resonance location (white dashed line 
in Fig.EIa)). 

In Fig. [9l^c) the 7/-4 resonance causes oscillations in 
ei that lead to escape when a > 0.71. In Fig. insta- 
bility is associated with the 5/-3 resonance separatrix. In 
Fig. [9l[e)&(f) the 3/-2 resonance causes oscillations in ei 
that lead to collision with the secondary when a > 0.835 
and a > 0.9, respectively. 

In Fig. [TO] we show the evolution of the 2/-1 resonant 
angle when /i = 0.05 for orbits with a = 0.6, ei = (a) 
and ei — 0.15 (b). When ei = the resonant angle circu- 
lates and when ei = 0.15 the resonant angle librates. The 
2/-1 resonant angle also exhibits short period oscillations 
which indicate that the assumption of Keplerian osculating 
elements in the short term is not very good even at mod- 
erate values of the mass ratio {/j, = 0.05). The surface of 
section (Fig llOf c)) confirms that these are regular orbits of 
the 2/-1 resonance, hence the spread of osculating elements 
is not due to chaotic diffusion. 



5.2 Retrograde case 

In Fig. [8] we show, for retrograde orbits, the maps with: (a) 
MEGNO chaos indicator; (b) times of disruption of 3-body 
system; (c) planet end states (stable, collision or escape). 

The location of the 2/-1 mean motion resonance, ob- 
tained by solving 5 = (Eq. ([9])) for a is shown as a white 
solid line in Figs. [SJa)&(b). For moderate to large fi values, 
Eq. (Illf) over-estimates the precession rate hence it displaces 
the theoretical resonance location to the left. We obtain a 
"corrected' 2/-1 mean motion resonance location by solving 
5 = for Q while taking into account the precession rate 
measured in the numerical integrations (this is shown as a 
white dashed line in Figs. [SJa)&(b)). 

The initial conditions that have C = Ci, C = C2 and 
C = are shown as black solid lines in Fig. [HI^c) (these 
curves approximately coincide). We know that the test par- 
ticle's end states depend on the the Jacobi constant. How- 
ever, although orbits in between the curves C = Ci, C = C2 
or C = can escape through Li, L2 or L3, respectively, 
in practice they only escape due to the effect of resonances. 
The stable region near a = 1 corresponds to test particles 



6 DISCUSSION 

We investigated the stability of prograde and retrograde 
planets in circular binary star systems. We saw that the 
cause of instability is either increase of eccentricity due to 
single mean motion resonance forcing, or chaotic diffusion of 
eccentricity and semi-major axis due to overlap of adjacent 
mean motion resonances. 

We computed the Jacobi constant in our grid of ini- 
tial conditions and compared with the values at Li, L2, 
I/3. We saw that the boundaries of the instability regions 
are explained by single resonance forcing or by resonances' 
overlap. Nevertheless, in the prograde planet's simulations, 
the ZVC opens at L\ near the instability border. However, 
in the retrograde planet's simulations, the ZVC opens at 
Li, L2 and L3 when a « 0.2 i.e. well before the instability 
border (a « 0.6). We conclude that the ZVC opening at 
Li, L2 or Lz are, as expected, necessary but not sufficient 
con ditions for instability . 

iQuarles et al.l l|201ll ) performed simulations of prograde 
orbits in binary systems and concluded that if the ZVC 
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Figure 9. Surfaces of section near stability boundary, (a) /i = 0.17: Orbits with a = 0.59 (red), a = 0.605 (magenta) and a = 0.61 
(blue) correspond, respectively, to inner circulation, separatrix, and outer circulation of the 2/1 retrograde resonance. The magenta orbit 
collides with the secondary after 28500 binary periods, (b) fi = 0.3: Orbit with a = 0.58 (red) exhibits inner circulation in the 2/1 
retrograde resonance. Orbit with a = 0.588 (blue) is at the separatrix of the 2/1 retrograde resonance and escapes after 600 binary 
periods, (c) /i = 0.175: Orbit with a = 0.7 (red) and a = 0.71 (blue) exhibit libration in 7/4 retrograde resonance, (d) fi = 0.14: Orbits 
with a = 0.74 (red) and a = 0.743 (blue) correspond, respectively, to inner circulation and separatrix of 5/3 retrograde resonance. The 
blue orbit escapes after 8890 binary periods, (e) /i = 0.08: Orbits with a = 0.8 (red), a = 0.82 (magenta) and a = 0.835 (blue) exhibit 
libration in 3/2 retrograde resonance, (f) fi = 0.05: Orbits with a = 0.89 (red) and a = 0.9 (blue) exhibit libration in 3/2 retrograde 
resonance. 



opens at L3 then the orbit is unstable. This does not contra- 
dict our results since our initial conditions differ. We start 
the planet at inferior conjuction as vi ewed from the c entra l 
star (i.e. between the 2 stars) while iQuarles et al.l (|201ll ) 
start the planet at opposition. In particular, in our prograde 
simulations the ZVC never opens at L3, while in our retro- 



grade simulations a large set of orbits with ZVC opening at 
L3 are stable. 

Prograde planets are unstable from a ~ 0.4 at /i > 0.15 
and from a ~ 0.5 at 0.15 > pi > 0.05. The main causes 
of instability are: overlap of the 4/1 resonance with nearby 
resonances at a ~ 0.4 and /j, > 0.15; overlap of the 3/1 reso- 
nance with nearby resonances at a ~ 0.5 and /i > 0.05; over- 
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Figure 8. Dynamical analysis for retrograde orbits following 
same criteria as in Fig. |4] a) Stability map in the («, fi) space 
showing < y >. b) Disruption times, c) End state of the planet. 
The theoretical 2/-1 resonance location shown in (a)&(b) as white 
solid line is over-displaced to the left when /i > 0.05. We show the 
"corrected" 2/-1 resonance location as white dashed line (see text 
for explanation). The black solid curves in (c) indicate the initial 
conditions with C = Ci, C = C2 and C = C3. The white squares 
at stability /instability transition zone indicate regions where we 
used the method of surfaces of section. The white circle is the 
initial condition for the red orbit in Fig. 1101 



Figure 10. Evolution of 2/1 resonance angle when fi = 0.05 and 
a = 0.6: ei = (a) and ei = 0.15 (b). Surface of section for both 
orbits (c). 

lap of the 5/2 and 2/1 resonances at a ~ 0.55 and ^ > 0.05; 
single resonance forcing in the 2/1 resonance at a ~ 0.6; 
overl ap of firs t order resonances when a > 0.7 in agreement 
with lWisdomI l)l980l ). 

Retrograde planets are unstable from a ~ 0.6 and fj, > 
0.15 which coincides with the location of the 2/-1 resonance. 
The instability border at a « 0.7 and fi — 0.175 coincides 
with the 7/-4 resonance. The instability border at a ~ 0.74 
and fj. — 0.14 coincides with the 5/-3 resonance separatrix. 
From a « 0.8 instability is due to eccentricity forcing at the 
3/-2 resonance^ 

[Mudryk fc Wul { 20061 ) identified the cause of prograde 
orbits' instability in eccentric binary systems as overlap of 
sub- resonances asso c iated with certain mean motion ratios 
p/q. iMudrvk fc Wul (|2006| ) also analyze previous low reso- 
lution numerical results of iHolman fc WiegertI l|l999t ) and 
argue that instability in circular binaries is due to overlap 
of sub-resonances associated with the 3/1 resonance. How- 
ever, this cannot work for circular binaries since there is 
only one resonant angle associated with a p/q resonance 
(all sub-resonances coincide). Here, we identified the cause 
of prograde orbits' instability in circular binary systems as 
single resonance forcing or overlap of different mean motion 
resonances, starting at the 4/1 resonance. 

We saw that retrograde planets are stable up to dis- 
tances closer to the perturber than prograde planets. We 
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conclude that this is due to essential difTerences between 
the phase-space topology of retrograde versus prograde res- 
onances. At mean motion ratio p/q, retrograde resonance 
has order p + q while prograde resonance has order p — q. 
Therefore, at a given resonance location a = {p/q)~^^^, we 
have: (a) eccentricity forcing on prograde planet is larger 
than eccentricity forcing on retrograde planet; (b) overlap 
with nearby resonances occurs at larger /i for retrograde 
con figuration than progra de configuration. 

iGavon fc Bo3 l|2008l ') showed, using MEGNO, that ret- 
rograde resonance in 2 planet systems is more stable than 
the equivalent prograde resonance. They conclude that this 
difference is due to close approaches being much faster and 
shorter for c ounter-revolving con figurations than for the pro- 
grade ones (|Gavon fc Boisll2008l ). While this is true, we be- 
lieve that the essential difference is not the duration of close 
approaches but instead, as explained above, the phase-space 
topology of prograde versus retrograde resonances. An ex- 
pansion of the Hamiltonian f or retrograde resonance in 2 
planet systems is presented in iGavon et al.l (|2009l ) but the 
numerical exploration of the model is limited to a small set 
of initial conditions and they do not conclude on the essen- 
tial differences between prograde and retrograde resonance. 

Finally, we refer that a similar mechanism could 
also explain the enhanced stability of retrograde satel- 
lites with respect to pr ograde satellite s that has been ob- 
served , for instance, by iHenonI (ll970i): [Hamilton fc Krivovl 
lll99'i1): iNesvornv et al.] l|2003h : Ishen fc Tremaind (|2008l ): 
iHinse et all (|2010l ). However, satellite motion is a distinct 
problem from that presented in this article (planet within 
binary system) as the hierarchy of masses is very different. 



APPENDIX A: 

Here, w e f ollow par tly the derivations i n iPeald (| 19761 ). 



WisdomI lll98Ct). iHenrard fc Lemaitrd USM) and 
Murrav fc DermottI (|l999l ') to model 1st. 2nd and 3rd 



order prograde resonances, and we extend this to model 3rd 
order retrograde resonance (2/-1). 

The Hamiltonian of the CR3BP near j/ (j — k) prograde 
or retrograde resonance is 



u__0__Hl,ri -LIT 

-n — -r fJres ~r *J sec 

2 ai 

where G mo — 1 ~ fi 
term is 



Ui af and G m2 



(Al) 

^. The resonant 



Ures = - — /d(a)ej cos((A; - j) Xi + j X2 - kvji) , (A2) 

0.2 



and the secular term is 



Us. 



a2 



(A3) 



From Lagrange's equations jMurrav fc DermottI [l999h . 
and Al change due to the secular term, while a\ and ei do 
not change; hence we can write 



H = — ^ „ h Ures 



2A2 

where we used Poincare canonical variables 

A = 111 , Al 

r = m ai(l - \/l - ej) , -Ti7i , 



(A4) 

(A5) 
(A6) 



At = (1 - ^1 - e2)ri7i* . 

Now, we change to resonant variables 

kK-{k-j)V 



= Al — A2 



and the secular variations in vj\ and Ai are 

^1 = 2 5 — y^') 

ni a\ e\ oe\ 

(A8) 
(A9) 

* = r ,1/. = [(fc-j)Ai+jA2-fcti7i]/fc (AlO) 
via the generating functioijfl 

F ^tp-H + (t)^ (All) 

Since F is time-dependent (A2 = ±712 t) the transformation 
introduces the term dF/dt in the new Hamiltonian, which 
is 



H 



, ^ -2 + Ures ± f 712 * T '^2 $ 



2 ($ 



(A12) 



Since H does not depend on 0, the momentum $ is a con- 
served quantity. Moreover, if ei ^ 1 then 'I' « "l>ef/2 and 
thus we expand the 1st term in H up to 2nd order around 
^ = 0. Hence, dropping constant terms that depend only on 
<I>, and changing the sign of H, we obtaiijf] 

H = 7* + ^*^ (-l)'''e(2*)''/^cos(A;V) (A13) 

wher^l 

7 = [{j -k){ni + Xl)Tjn2 + kTbl]/k (A14) 
3 (j - kf 



2 k^af 



r / \ —k/2 —k 

e = — }d{OL)n-^ ai 

0.2 



(A15) 
(A16) 



As explained in lMurrav fc DermottI (|l999l ). by introducing 
a scaled momentum 



(A17) 



this can be written as a single parameter Hamiltonian 

H ^ 5^ + ^!'^ + 2(~lf {2^)''''^ cos{kil}) (A18) 

with 



5 = 7 



£2 ^2 



(A19) 



The Hamiltonian (Eq. IA18|l can be expressed in cartesian 
canonical variables {x — R cos{ip),y = R smitp)) where the 
scaling factor is 



R 



2^ . 



(A20) 



For a mixed variable transformation: (p = dF/d^, V = dF/d^, 
A = dF/dXi, r = -dF/d'uji. 

* Where (—1)''" is introduced because e > if fc even and e < if 
k odd. 

^ Not e that /3 differs from expression in iMurrav &: Dermotd 
1I1999I) (no mass). 
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